#  R3.6

reference:Controlled modelling of human epiblast and amnion development using stem cells

# loading R library
#R3.6 
rm(list=ls())
rewrite=FALSE

condaENV <- "/home/chenzh/miniconda3/envs/R3.6"
LBpath <- paste0(condaENV ,"/lib/R/library")
.libPaths(LBpath)

suppressMessages(library(dplyr))
suppressMessages(library(tidyr))
suppressMessages(library(ggplot2))
suppressMessages(library(cowplot))
suppressMessages(library(pheatmap))
suppressMessages(library(Seurat))



# working directory
DIR <- "~/My_project/CheckBlastoids"
setwd(DIR)
knitr::opts_knit$set(root.dir=DIR)

Loading R functions

#source("~/PC/R_code/functions.R")
#source("~/PC/SnkM/SgCell.R")
source("src/local.quick.fun.R")
feature.plot.cor <- c("yellow","red","black")

suppressMessages(library(foreach))
suppressMessages(library(doParallel))
numCores <- 10
registerDoParallel(numCores)

options(digits = 4)
options(future.globals.maxSize= 3001289600)
TD="Nov14_2021"

loading gene metadata

load("tmp_data/gene.meta.Rdata",verbose=T)
## Loading objects:
##   gtf.anno
##   mt.gene
##   ribo.gene
savefile=paste0("tmp_data/",TD,"/EPSC_Amnion.Rdata")

if (file.exists(savefile)) {
  load(savefile,verbose = T)
}else{
  load(paste0("tmp_data/",TD,"/all.counts.meta.Rdata"),verbose=T)
  
 
  #' loading published seurat object
  AMLC.ob <- readRDS("data/GSE134571/GSE134571_Posterior48h_H9_Amnion_Merged.rds")
  AMLC.ob <- AMLC.ob %>% RenameIdents('0'="Tsw-AMLC",'1'='MeLC2','2'='hES','3'='MeLC1','4'='hPGCLC','5'='AMLC')
  
  AMLC.DM <- FunRF_FindAllMarkers_para(AMLC.ob)
  temp.plot1 <- list()
  temp.plot1$id <- DimPlot(AMLC.ob,label=T) +NoAxes()+NoLegend()
  temp.plot1$oid <- DimPlot(AMLC.ob,label=T,group.by="orig.ident") +NoAxes()+NoLegend()
  temp.mk.gene <- (AMLC.DM$sig  %>% filter(NP=="pos") %>% filter(gene %in% rownames(counts.all)) %>% group_by(set) %>% top_n(2,power) %>% pull(gene))
  for (g in temp.mk.gene) {
    temp.plot1[[g]] <- FeaturePlot(AMLC.ob,g)+NoAxes()+NoLegend()+ggtitle(g)+FunTitle()
  }
  
  #same QC from the paper
  temp.raw <-  meta.all %>% filter(pj=="JPF2019")
  temp.M <- meta.all %>% filter(pj=="JPF2019" & EML=="Pos48ELS") %>% filter(nGene >=3200 & nGene <= 6200 & mt.perc < 0.06) %>% bind_rows(meta.all %>% filter(pj=="JPF2019" & EML=="H9AMN") %>% filter(nGene >=3600 & nGene <= 6400 & mt.perc < 0.06))
 
  counts.filter <- counts.all[setdiff(rownames(counts.all), mt.gene),temp.M$cell]
  DR.TPM.filter <- apply(counts.filter,2,function(x){x/sum(x)*1000000})
  
  # expressed genes
  temp.sel.expG <- rownames(counts.filter)[rowSums(counts.filter[,c(temp.M$cell)] >=1) >5]
  
  data.temp <-  CreateSeuratObject(counts.filter[temp.sel.expG,temp.M$cell], meta.data = (temp.M %>% tibble::column_to_rownames("cell"))) %>% NormalizeData(verbose = FALSE) %>% FindVariableFeatures( selection.method = "vst", nfeatures = 2000, verbose = FALSE)  %>% ScaleData(verbose=F,vars.to.regress=c("mt.perc"))%>% RunPCA(verbose=F,npcs=20) %>% RunUMAP(dims=1:20,verbose=F,min.dist=0.4) %>% FindNeighbors( dims = 1:20,verbose = FALSE) %>%  FindClusters(resolution = 0.45,verbose = FALSE)  #"mt.perc",
  data.mod <- data.temp
  id <- Idents(data.mod) %>% as.vector()
  names(id) <- colnames(data.mod)
  id[colnames(data.mod)[data.mod@meta.data$EML=="H9AMN" & (id %in% c(0,5))]] <- "Tsw-AMLC"
  id[colnames(data.mod)[data.mod@meta.data$EML=="H9AMN" & (id %in% c(1))]] <- "hES"
  id[colnames(data.mod)[data.mod@meta.data$EML=="Pos48ELS" & (id %in% c(3))]] <- "hPGCLC"
  id[colnames(data.mod)[data.mod@meta.data$EML=="Pos48ELS" & (id %in% c(2))]] <- "MeLC1"
  id[colnames(data.mod)[data.mod@meta.data$EML=="Pos48ELS" & (id %in% c(4,7))]] <- "MeLC2"
  id[colnames(data.mod)[data.mod@meta.data$EML=="Pos48ELS" & (id %in% c(6))]] <- "AMLC"
  id[! id %in% c("Tsw-AMLC","hES","hPGCLC","MeLC1","MeLC2","AMLC")] <- "mix"
  Idents(data.mod) <- as.factor(id)
  data.mod <- subset(data.mod,cells=names(id)[id!="mix"])
  
  table(data.mod %>% Idents())

  save(data.mod,temp.sel.expG,temp.M,temp.raw,temp.plot1,AMLC.DM,temp.mk.gene,file=savefile)
}
## Loading objects:
##   data.mod
##   temp.sel.expG
##   temp.M
##   temp.raw
##   temp.plot1
##   AMLC.DM
##   temp.mk.gene

repeat paper’s figure

cowplot::plot_grid(plotlist=temp.plot1)

my new annotation

temp.plot2 <- list()
temp.plot2$id <- DimPlot(data.mod,label=T) +NoAxes()+NoLegend()
## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.
temp.plot2$oid <- DimPlot(data.mod,label=T,group.by="EML") +NoAxes()+NoLegend()
for (g in temp.mk.gene) {
  temp.plot2[[g]] <- FeaturePlot(data.mod,g)+NoAxes()+NoLegend()+ggtitle(g)+FunTitle()
}
cowplot::plot_grid(plotlist=temp.plot2)

#VlnPlot(data.mod,temp.mk.gene)

General distribution of nGene and mt.perc

temp.plot <- list()
temp.plot$AMN1 <- temp.raw %>% filter(EML=="Pos48ELS")%>% ggplot +geom_point(mapping=aes(x=mt.perc,y=nGene,col=EML))+geom_vline(xintercept = 0.06,linetype = "dashed" )+geom_hline(yintercept = 6200,linetype = "dashed")+geom_hline(yintercept = 3200,linetype = "dashed")+theme_classic()+ggtitle("QC for 10X(Pos48ELS)")+FunTitle()
temp.plot$AMN2 <- temp.raw %>% filter(EML=="H9AMN")%>% ggplot +geom_point(mapping=aes(x=mt.perc,y=nGene,col=EML))+geom_vline(xintercept = 0.06,linetype = "dashed" )+geom_hline(yintercept = 6400,linetype = "dashed")+geom_hline(yintercept = 3600,linetype = "dashed")+theme_classic()+ggtitle("QC for 10X(H9AMN)")+FunTitle()
cowplot::plot_grid(plotlist=temp.plot,nrow=2,ncol=2)

saving the new meta data

AMN.meta <- temp.M %>% filter(cell %in% colnames(data.mod)) 
AMN.meta$EML <- Idents(data.mod)[AMN.meta$cell] %>% as.vector()

downsample dataset

FunMaSF <- function(temp,n) {
  set.seed(123)
  return(head(temp[sample(nrow(temp)),],n))
}
temp.out <- NULL

meta.AMN.further.ds <- split(AMN.meta,AMN.meta$EML) %>% lapply(function(x){return(FunMaSF(x,100))}) %>% do.call("bind_rows",.)

if (!file.exists(paste0("tmp_data/",TD,"/meta.AMN.further.rds")) | rewrite) {
  print("save output")
  saveRDS(AMN.meta,paste0("tmp_data/",TD,"/meta.AMN.further.rds"))
  saveRDS(meta.AMN.further.ds,file=paste0("tmp_data/",TD,"/meta.AMN.further.ds.rds"))
}

sessionInfo()                    
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server 7.9 (Maipo)
## 
## Matrix products: default
## BLAS/LAPACK: /home/chenzh/miniconda3/envs/R3.6/lib/libopenblasp-r0.3.7.so
## 
## locale:
##  [1] LC_CTYPE=en_US.utf-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.utf-8        LC_COLLATE=en_US.utf-8    
##  [5] LC_MONETARY=en_US.utf-8    LC_MESSAGES=en_US.utf-8   
##  [7] LC_PAPER=en_US.utf-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.utf-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] doParallel_1.0.15 iterators_1.0.12  foreach_1.4.8     Seurat_3.1.4     
## [5] pheatmap_1.0.12   cowplot_1.0.0     ggplot2_3.2.1     tidyr_1.1.0      
## [9] dplyr_1.0.0      
## 
## loaded via a namespace (and not attached):
##   [1] TH.data_1.0-10      Rtsne_0.15          colorspace_1.4-1   
##   [4] ellipsis_0.3.1      ggridges_0.5.2      farver_2.0.3       
##   [7] leiden_0.3.3        listenv_0.8.0       npsurv_0.4-0       
##  [10] ggrepel_0.8.2       mvtnorm_1.0-12      codetools_0.2-16   
##  [13] splines_3.6.2       mnormt_1.5-6        lsei_1.2-0         
##  [16] knitr_1.22          TFisher_0.2.0       jsonlite_1.6.1     
##  [19] ica_1.0-2           cluster_2.0.8       png_0.1-7          
##  [22] uwot_0.1.5          sctransform_0.2.1   compiler_3.6.2     
##  [25] httr_1.4.1          Matrix_1.2-17       lazyeval_0.2.2     
##  [28] htmltools_0.3.6     tools_3.6.2         rsvd_1.0.3         
##  [31] igraph_1.2.5        gtable_0.3.0        glue_1.4.1         
##  [34] RANN_2.6.1          reshape2_1.4.4      rappdirs_0.3.1     
##  [37] Rcpp_1.0.4.6        Biobase_2.46.0      vctrs_0.3.0        
##  [40] multtest_2.42.0     gdata_2.18.0        ape_5.3            
##  [43] nlme_3.1-139        gbRd_0.4-11         lmtest_0.9-37      
##  [46] xfun_0.6            stringr_1.4.0       globals_0.12.5     
##  [49] lifecycle_0.2.0     irlba_2.3.3         gtools_3.8.2       
##  [52] future_1.16.0       MASS_7.3-51.3       zoo_1.8-8          
##  [55] scales_1.1.1        sandwich_2.5-1      RColorBrewer_1.1-2 
##  [58] yaml_2.2.0          reticulate_1.14     pbapply_1.4-2      
##  [61] gridExtra_2.3       stringi_1.4.6       highr_0.8          
##  [64] mutoss_0.1-12       plotrix_3.7-8       caTools_1.18.0     
##  [67] BiocGenerics_0.32.0 bibtex_0.4.2.2      Rdpack_0.11-1      
##  [70] rlang_0.4.6         pkgconfig_2.0.3     bitops_1.0-6       
##  [73] evaluate_0.13       lattice_0.20-38     ROCR_1.0-7         
##  [76] purrr_0.3.4         labeling_0.3        patchwork_1.0.0    
##  [79] htmlwidgets_1.3     tidyselect_1.1.0    RcppAnnoy_0.0.15   
##  [82] plyr_1.8.6          magrittr_1.5        R6_2.4.1           
##  [85] gplots_3.0.3        generics_0.0.2      multcomp_1.4-12    
##  [88] pillar_1.4.4        withr_2.2.0         sn_1.5-5           
##  [91] fitdistrplus_1.0-14 survival_3.1-12     tibble_3.0.1       
##  [94] future.apply_1.4.0  tsne_0.1-3          crayon_1.3.4       
##  [97] KernSmooth_2.23-15  plotly_4.9.2        rmarkdown_2.3      
## [100] grid_3.6.2          data.table_1.14.0   metap_1.3          
## [103] digest_0.6.25       numDeriv_2016.8-1.1 RcppParallel_4.4.4 
## [106] stats4_3.6.2        munsell_0.5.0       viridisLite_0.3.0